Electronic compressibility of a graphene bilayer 



S. Viola Kusminskiy, Johan Nilsson, D. K. Campbell, and A. H. Castro Neto 1 
1 Department of Physics, Boston University, 590 Commonwealth Ave., Boston, MA 02215 

We calculate the electronic compressibility arising from electron-electron interactions for a 
graphene bilayer within the Hartree-Fock approximation. We show that, due to the chiral na- 
ture of the particles in this system, the compressibility is rather different from those of either the 
two-dimensional electron gas or ordinary semiconductors. We find that an inherent competition be- 
tween the contributions coming from intra-band exchange interactions (dominant at low densities) 
and inter-band interactions (dominant at moderate densities) leads to a non-monotonic behavior of 
the compressibility as a function of carrier density. 
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The recently developed experimental capability of iso- 
lating and manipulating an arbitrary number of graphene 
layers [l[ has attracted considerable attention both for 
its impact on basic science Q and for the tantalizing po- 
tential technological applications. The graphene bilayer 
is particularly interesting because of the possibility of 
opening - and controlling - a gap in the electronic spec- 
trum by applying an external electric field 0, 0, [a @] • 
This is not possible for the single layer graphene. The 
bilayer, therefore, while inheriting many of the peculiar 
electronic characteristics of the monolayer due to its chi- 
ral Dirac fermion (though massive) spectrum, has the 
added virtue of being capable of acting as an electronic 
switch. It is thus essential to obtain a comprehensive 
characterization of this material. While some transport 
experiments are available Q, thermodynamic measure- 
ments are largely lacking. Among the thermodynamic 
quantities to be measured, the electronic compressibility 
k stands out as an excellent tool to provide insight into 
the many-body interactions present in this material, n 
can be obtained from the ground state energy as: 



k- 1 = n 2 e {d 2 E/dn 2 e ). 



(1) 



where E is the ground state energy per unit area, and n e 
is the electronic density. The electronic compressibility 
of a single layer graphene has been recently measured [8( , 
and its behavior, besides being remarkably different from 
that of the usual two-dimensional gas (2DEG), seems to 
indicate that contributions from Coulomb interactions 
are either very weak or cancel out. Hartree-Fock 0] and 
Random Phase Approximation (RPA) [lOJ] calculations 
predict a correction between 10% and 20% to the free 
theory for experimentally realized dopings. This correc- 
tion increases logarithmically as the doping is lowered. It 
is natural then to ask what role interactions play in the 
bilayer. In many aspects, the bilayer graphene closely 
resembles the 2DEG, as described below. Hence, the bi- 
layer system provides an opportunity to isolate the effects 
arising from its single layer constituents, from those oc- 
curring in an ordinary 2DEG. In particular, the issue of 
the chirality, which is so important for weak-localization 
physics [ll|], is the main difference between these two 



systems, and as we will show, plays an important role 
in the many-body physics of the bilayer. For small dop- 
ing, the bilayer can be mapped approximately to a chiral 
two-dimensional massive fermionic system with parabolic 
bands [H, EH, in contrast with the massless, cone-like 
dispersion found in the monolayer. This limit is useful to 
compare the behavior of the bilayer (with its chirality) 
to that of the ordinary 2DEG, where experiments have 
shown that interactions play a dominant role, making the 
proper compressibility negative for small electron densi- 
ties [14( as opposed to a positive constant given by the 
non-interacting model. This behavior is already present 
at the Hartree-Fock level. Due to the aforementioned 
mapping, a priori it is reasonable to expect that the ef- 
fect of electron-electron interactions would be observable 
in the graphene bilayer. 

In this paper we calculate within the Hartree-Fock ap- 
proximation the dependence of the inverse compressibil- 
ity on the Fermi vector kp using both a full, four-band 
(4B) model and the two-band (2B) approximation, which 
is valid for very small doping. We show that the most im- 
portant qualitative signatures of the compressibility are 
already present at the 2B model level but that the 4B cal- 
culation, while more cumbersome, reveals finer features. 

Throughout this paper we will refer (loosely) to the 
quantity k~ l = dn/dn e as the "inverse compressibility". 
Here /x stands for the chemical potential of the system. 
k differs by a factor of n 2 from k in {2}. This is ap- 
propriate since k^ 1 is usually the actual experimentally 
measured quantity. The density of electrons is given by 
n e = gsg v kp/ (^ir), with g s — 2, g v = 2 being the spin 
and valley degeneracy, respectively In the following we 
will consider the case of small doping but outside the 
range of ferromagnetic instability that is found at ex- 
tremely low doping [13]. At the Hartree-Fock level, the 
ground state energy is given by E = K + E ex , where 
K stands for the kinetic energy and E ex is the exchange 
energy per unit of area. 

A graphene bilayer consists of two planes of graphene 
stacked as shown in Fig. (p}. The kinetic term 
of the Hamiltonian can be written, in the near- 
est neighbor tight binding approximation [l5| by 
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in the continuum as: 



FIG. 1: (Color online) Graphene bilayer lattice showing the 
two underlying sublattices A (purple) and B (yellow). In- 
set: 4B dispersion (continuous line) and 2B approximation 
(dashed) as a function of momentum k. Momentum is in 
units of the cutoff A, measured from the K point. 



expanding around the K, K' points of the Bril- 
louin zone, as TLkin = J2q V>q£(p)^Q witn ^>q = 



( c l a ^„, c 1r ~„,cl R „ „ , c„ a „ „ ) where a labels the 
valley, a the spin and Ai, Bi denotes the sublattice in 
the plane i = 1,2. J2q represents the sum over all the 
indices. The kinetic energy matrix is given by (we use 
units such that h = 1): 
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where tan</>(p) = p y /p x , t± = 0.35 eV is the inter-layer 
hopping energy and we have set v F = 3to/2 6.6 eV A) 
to unity (t is the intra-layer hopping energy and a the in- 
plane carbon-carbon distance). The interaction is given 
by the 2D Fourier transform of the 3D Coulomb poten- 
tial, which is Vi P (k) = ^p-^ for the interaction among 

electrons within the same plane and V^k) = 2 ^r—^- k — 
otherwise, being d ~ 3.35 A the inter-plane distance. 

The kinetic energy matrix ([2]) can be diagonalized by a 
unitary transformation <S^(p). The resulting dispersion 
bands (see Fig.Q) are: Ei(p) = -t + E(p), E 2 (p) = 
t-E(p), E 3 (p) = t + E(p) and E 4 (p) = -t-E(p); being 
E(p) — \/t 2 + p 2 and t — t±/2. It is convenient to work 
with the symmetric and anti-symmetric combinations of 
the layer densities, p± — pi±p 2 , which can be expressed 
in the diagonal basis as p a (<l) = X^p^^P + q)x"(P + 
q, p^p) with $(p) = S^(p)i;(p) and a = ±. The 4x4 
matrices x a contain the information of the overlap due 
to the change of basis. Then the interaction Hamiltonian 
takes the form Hi = 1/(2 A) J2q^o J2 a Pa(q)K*(qW(q), 
and the exchange energy per unit area A can be written 



1 f d 2 p ePq ^ 

E ex = -g s g v - J (2^)2 (2^)2 L ^(*P)^(P.q 



ni(q)nj(p)V a (q- p) , 



(3) 



where i, j = 1, ...4; n t (q) =< $t( q )$ i ( q ) > and V±(k) = 
i(Vi P (k) ± V op (k)). The occupation factors are given 
by nx(q) = <d(k F - q) [n x (q) = 0], n 2 (q) = 1 [n 2 (q) 
I — <d(kp — q)], n 3 (q) = and 714(g) = 1 in the case of 
electron [hole] doping. This model however requires a 
cutoff A of the order of the inverse of the lattice param- 
eter. 

Being simpler to work with, and widely used as a start- 
ing point for calculations in the graphene bilayer, we start 
our analysis with the approximate 2B model that can 
be constructed at low energies by performing degenerate 
perturbation theory 1^, HI]. This results in an effective 
kinetic Hamiltonian: 
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with = (c^ Bl .a^ c t,A 2 ^,a)- The result of the ap- 
proximation is an effective model with opposite parabolic 
dispersion bands of energy E a = p 2 /(2t), Eh — —p 2 /(2t) 
as shown in Fig. (JT|) . The effective kinetic energy per 
unit area then is given by K = (k F — A 4 )/(47rf) giv- 
ing a kinetic contribution to the inverse compressibility 
n/(2t). 

In this reduced Hilbert space Eq. ([3]) is still valid, but 
this time the x"(Piq) are 2x2 matrices. Combining all 
the contributions and re-inserting the units, we find the 
total inverse compressibility k _1 in the 2B model to be 
given by the expression: 
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Here we have defined r(p, q 7 9) — \/p 2 + q 2 — 2pq cos 9, t, 
k F and 1/d are in units of A, and g = e 2 /(v F eo) is the 
graphene coupling strength. The ± indicates the expres- 
sion for electrons (+) or holes (— ). The differing term 
however is roughly a constant and can be neglected for 
small doping. Therefore, in what follows we will use the 
results for electron doping. 

Fig. ((21) shows a plot of ([5| as a function of k F for 
i/A = 0.026, and dA = 3.7 (A = 1.06 w 1 A" 1 ). As can 
be seen, for very small doping the compressibility changes 
sign, becoming negative and divergent. This behavior, as 
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FIG. 2: (Color online) Inverse compressibility vs. Fermi wave 
vector in the 2-band approximation. Total: solid line (red), 
intra-band plus kinetic contribution: dashed (blue), kinetic: 
dash/dot (green), 2DEG: dotted (black). The inset depicts 
the contribution of intra-band exchange (solid line) and inter- 
band exchange (dashed). 



mentioned previously, is also observed in the 2DEG. It is 
instructive to discriminate between inter and intra-band 
contributions. Fig.@ also depicts the inverse compress- 
ibility when only the intra-band transitions are consid- 
ered (as well as the kinetic term, of course). From the 
difference with the curve for the total inverse compress- 
ibility, we can conclude that the inter-band contribution 
tends to move the negative region to smaller densities. 
The overall effect of the inter-band transitions then is 
to enhance therefore reducing the compressibility. 
This can be seen clearly from the inset in Fig.@, where 
the contribution from intra-band transitions is negative 
while that from inter-band transitions is positive. Apart 
from the sign, both present a similar behavior and are 
comparable in magnitude. In the 2B approximation, the 
kinetic contribution, as in the 2DEG case, is independent 
of the electronic density and is also plotted in Fig. @ for 
reference. Therefore, the resulting total compressibility 
will be given by a competition of the two contributions. 
This difference in sign between inter and intra-band con- 
tributions is analogous to the one present in monolayer 
graphene [8| and, whitin the 2B approximation, it can 
be interpreted in terms of the chirality of the quasiparti- 
cles. Intra (inter)-band exchange corresponds to interac- 
tions between particles of the same (opposite) chirality. 
Remarkably though, while for the monolayer the total 
exchange contribution is positive, for the bilayer is nega- 
tive. 

We can also compare with the usual result for the 
2DEG. For this we start from the expression for the chem- 
ical potential (see for example [161]) fi = d[ne(n)]/dn, 
where n = k F /ir is the electronic density (taking into 
account spin and valley degeneracy) and e is the ground 
state energy per electron. If we consider only kinetic 



and exchange energy, e = eq + e ex and for the 2DEG 
e o = kp/(4m), e ex = —4e 2 k F Therefore k^deg = 
(n/m — 2e 2 / kp)/2. To compare with our case, we identify 
t = m and write k-^deg ~ u -f/(2A) ~ 2<7/&f). The 
corresponding plot is also depicted in Fig. @ . Within the 
2B approximation, the bilayer compressibility behaves 
qualitatively similar to that of the 2DEG, although, be- 
cause of the chiral nature of the bilayer system, the region 
of negative values is shifted to smaller values of the Fermi 
vector and therefore smaller densities. 

As mentioned above, the total compressibility of the 
bilayer is a product of the competition between inter 
and intra-band contributions to exchange. However from 
the inset in Fig. |(3J) it transpires that this is a very tight 
competition and small changes may alter the result. We 
therefore proceed to analyze the full 4B model. Due to 
the symmetries of the x matrices, expression ([3J is greatly 
simplified. The exchange contribution to the compress- 
ibility then is given by 
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where 9 is the angle between kp, p. As the notation sug- 
gests, -D++ corresponds to exchange within the positive 
conduction band 1 while D? measures the exchange 
between the negative filled sea and the conduction band. 
The calculation for hole doping is completely analogous, 
with the overlap elements to be considered for that case 
being jx^l^F, p)| 2 and |x24(kF, p)| 2 - Since the com- 
pressibility involves only occupied states, its behavior 
is not symmetric with respect to particle-hole exchange. 
Nonetheless, the explicit calculation shows that the dif- 
ference is negligible for small doping, being the same as 
in the 2B case. On the other hand, the kinetic con- 
tribution to the inverse compressibility is independent 
of the type of carrier and it is easily calculated to be 
k^ 1 = n/[2E(p)]. The final results for the four band cal- 
culation, by summing all the contributions, is shown in 
Fig. ©. 

We see that the full model confirms the major quali- 
tative features found within the 2B approximation. The 
compressibility is negative for small electronic density, 
diverging in the limit of k F — > 0, and the inter-band ex- 
change contributes to the incompressibility of the system 
(see Fig. flU). The difference observed at larger values 
of the Fermi momentum is simply due to the difference 
in the kinetic term, since it is a constant in the two band 
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FIG. 3: (Color online) Inverse compressibility as a function 
of the Fermi wave vector as calculated from the 4B model. 
Negative values of the Fermi vector indicate the result for 
hole doping. 
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FIG. 4: (Color online) Comparison between the exchange con- 
tribution for the 2B and 4B models. Total exchange (red): 4B 
solid line, 2B dotted; intra-band only (blue): 4B dashed, 2B 
dash/dot. 



portant to explain the behavior of the compressibility. 
The non-monotonic behavior of the compressibility ob- 
tained with the 4B model is highly unusual. Generally, 
nonmonotonicity is associated with some external factor, 
such as confinement or applied magnetic fields. However 
here it is due solely to intrinsic electronic interactions. 
The implications of this remain to be understood. On 
the other hand, the negativity of the compressibility is 
understood once it is realized that the one involved is 
not the total compressibility but only that of the elec- 
tronic gas. The total compressibility will comprise also 
the positive ionic background, which stabilizes the sys- 
tem. Nonetheless, the negative divergence of the inverse 
compressibility, present in both models for low enough 
electronic densities, could signal the eventual onset of 
Wigner crystallization [T^. These results, as in the case 
of the single layer graphene [|| can be studied via sin- 
gle electron transistor (SET) measurements. Our results 
indicate that the compressibility turns negative at den- 
sity values of approximately n e rj 10 11 /cm 2 , which bor- 
ders the current available precision However, being a 
Hartree-Fock calculation, this can act just as a very rough 
estimate. We have also neglected the trigonal warping 
term, which might be of importance at very low densi- 
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case while in the 4B model it is ~ I/^f for large hp- 
However, a more detailed comparison reveals a peak in 
the inverse compressibility that is not captured in the 2B 
model. Overall, the 4B calculation predicts a more in- 
compressible system in the range for which k > 0, up to 
approximately three times larger than the prediction of 
the 2B model at the position of the peak. As seen from 
Fig. Q , the bulk of the difference between the 2B result 
and the 4B one comes from the inter-band exchange. 

The 4B model seems to predict a behavior that is a hy- 
brid between the one of a 2DEG, where the total contri- 
bution to the compressibility from exchange is negative, 
and the graphene monolayer, where it is positive. 

In conclusion, we have studied the electronic compress- 
ibility of a graphene bilaycr within the Hartree-Fock ap- 
proximation and have found a behavior that is remark- 
ably different from the two-dimensional electron gas due 
to the presence of inter-band transitions, and also from 
graphene monolayer. We have shown that the inverse 
compressibility is not a monotonic function of the elec- 
tronic density and that the effective 2B model gives a 
good description of the problem only at very low densi- 
ties. At intermediate densities, the four bands are im- 
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